# Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
#### 1. 课堂演示-衍生得到相应的变量，同步衍生 ####
##### 1.1 性别-Sex#####
paper.data$Sex <- ifelse(paper.data$RIAGENDR == 1, 'male', 'female')
table(paper.data$Sex)

tab1 <- tableby(Type ~ factor(Sex),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)
##### 1.2 年龄-age.group #####
paper.data$age.group <- ifelse(paper.data$RIDAGEYR >= 60 & paper.data$RIDAGEYR < 69, '60-69 years',
                               ifelse(paper.data$RIDAGEYR >=70 & paper.data$RIDAGEYR < 79, '70-79 years',
                                      '80+ years'))
table(paper.data$age.group)

tab1 <- tableby(Type ~ factor(age.group),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)

##### 1.3 种族-race #####
table(paper.data$RIDRETH3)
race <- recode_factor(paper.data$RIDRETH3, 
                      `1` = 'Mexican American',
                      `2` = 'Other Hispanic',
                      `3` = 'Non-Hispanic White',
                      `4` = 'Non-Hispanic Black',
                      `6` = 'Other/multiracial',
                      `7` = 'Other/multiracial'
                      )
paper.data$race <- race
table(paper.data$race)

tab1 <- tableby(Type ~ factor(race),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)


##### 1.4 复杂转换-饮酒 alq.group #####
# 仅衍生有无饮酒
# 区分男性和女性不同的饮酒标准，drinks per day (dkspd) ：
#### 是否酗酒：https://bmcgastroenterol.biomedcentral.com/articles/10.1186/s12876-022-02430-7
#### 涉及男性和女性的饮酒标准不一致：https://www.nature.com/articles/s41598-017-02426-4  


# 衍生方式：Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# 衍生结果为：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/ALQ_G.htm#ALQ120U （课程使用的paper）
# ALQ120Q, 具体的数值;  777:refused, 999:don't know
# ALQ120U, 1:week, 2:month, 3:year, 7:refused, 9:don't know

# 使用 plyr 包的 ddply 快速看一下数据
ddply(paper.data, .(ALQ101, ALQ110), summarize, n = length(SEQN))

ddply(alq.g, .(ALQ101, ALQ110), summarize, n = length(SEQN))

# 统一单位为月
# week-month: *4; year-month:/12
ori.alq.unit <- paper.data$ALQ120U
table(ori.alq.unit)
trans.unit.month <- ifelse(ori.alq.unit == 1, 4, 
                           ifelse(ori.alq.unit == 3, 1/12, 
                                  ifelse(ori.alq.unit == 7|ori.alq.unit == 9, NA, 1)))
paper.data$trans.unit.month <- trans.unit.month
# View(paper.data)


# 统一数值为月饮酒量-quantity
ori.alq.quantity <- paper.data$ALQ120Q
trans.quantity.month <- ifelse(ori.alq.quantity >= 0,  
                               ori.alq.quantity * trans.unit.month, NA)
paper.data$trans.quantity.month <- trans.quantity.month
View(paper.data[,c('SEQN','ALQ101','ALQ110', 'ALQ120Q', 'ALQ120U', 'trans.quantity.month', 'trans.unit.month')])


# 将饮酒量的结果划分为4档：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# 首先根据转化得到的trans.quantity.month，将饮酒量划分3档
# 利用 paper.data$ALQ101
alq101 <- paper.data$ALQ101
trans.quantity.month.factor <- ifelse(trans.quantity.month >=1 & trans.quantity.month <5, '1-5 drinks/month',
                                      ifelse(trans.quantity.month >=5 & trans.quantity.month <10, '5-10 drinks/month',
                                             ifelse(trans.quantity.month >=10, '10+ drinks/month', 'wait')))

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor

ddply(paper.data, .(ALQ101, trans.quantity.month.factor), summarise, n = length(SEQN))
#View(paper.data[,c('SEQN','ALQ101','ALQ110', 'ALQ120Q', 'ALQ120U', 'trans.quantity.month', 'trans.unit.month', 'trans.quantity.month.factor')])

# ALQ101 回答是1，但是 trans.quantity.month.factor 没有具体的值，则将 trans.quantity.month.factor 取值为1
table(trans.quantity.month.factor) 
index.1 <- which((trans.quantity.month.factor=='wait' | is.na(trans.quantity.month.factor)) & alq101 == 1) #844
trans.quantity.month.factor[index.1] <- '1-5 drinks/month'
table(trans.quantity.month.factor)

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor
ddply(paper.data, .(ALQ101, trans.quantity.month.factor), summarise, n = length(SEQN))

index.less.1 <- which(alq101 == 2)
trans.quantity.month.factor[index.less.1] <- 'Non-drinker'

table(trans.quantity.month.factor)
# 根据是否饮酒补充 Non-drinker

paper.data$alq.group <- trans.quantity.month.factor
# View(paper.data)

table(paper.data$alq.group)

tab1 <- tableby(Type ~ factor(alq.group) + factor(RIAGENDR),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)
# View(paper.data[,c('SEQN','ALQ101','ALQ110', 'ALQ120Q', 'ALQ120U',
#                      'trans.unit.month', 'trans.quantity.month','alq.group')])



##### 1.5 Total 卡路里-total.calories #####
# day1、day2 都有，则取2天平均，否则，取第1天的值
# View(paper.data[,c('DR1TKCAL', 'DR2TKCAL')])
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.calories <- apply(paper.data[,c('DR1TKCAL', 'DR2TKCAL')], 1, mean)
paper.data$total.calories <- total.calories
#View(paper.data[,c('DR1TKCAL', 'DR2TKCAL', 'total.calories')])

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$DR2TKCAL))
total.calories[day.2.na.index] <- paper.data$DR1TKCAL[day.2.na.index]
total.calories[day.2.na.index]
#View(paper.data[,c('DR1TKCAL', 'DR2TKCAL', 'total.calories')])

# 添加到paper.data中，再次确认
paper.data$total.calories <- total.calories
# View(paper.data[,c('DR1TKCAL', 'DR2TKCAL', 'total.calories')])


##### 1.6 Total 饮食中 L&Z-total.dr.LZ#####
# https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/DR2TOT_G.htm#DR2TLZ
# day1、day2 都有，则取2天平均，否则，取第一天的值
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.dr.LZ <- apply(paper.data[,c('DR1TLZ', 'DR2TLZ')], 1, mean)
paper.data$total.dr.LZ <- total.dr.LZ
#View(paper.data[,c('DR1TLZ', 'DR2TLZ', 'total.dr.LZ')])

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$total.dr.LZ))
total.dr.LZ[day.2.na.index] <- paper.data$DR1TLZ[day.2.na.index]
total.dr.LZ[day.2.na.index]

# 添加到paper.data中，再次确认
paper.data$total.dr.LZ <- total.dr.LZ
# View(paper.data[,c('DR1TLZ', 'DR2TLZ', 'total.dr.LZ')])

##### 1.7 认知评分 #####
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3

colnames(paper.data)
CERAD.total <- apply(paper.data[,c('CFDCST1', 'CFDCST2', 'CFDCST3')], 1, sum)
paper.data$CERAD.total <- CERAD.total
# View(paper.data[,c('CFDCST1', 'CFDCST2', 'CFDCST3','CERAD.total')])

#### 2. 变量Summary ####
nhanes.design <- svydesign(data=paper.data[which(paper.data$WTDRD1 > 0),], 
                           id=~SDMVPSU, strata=~SDMVSTRA, 
                           weights=~WTDRD1, nest=TRUE)
summary(nhanes.design)

##### 2.1 离散变量 Summary #####
svytotal(~age.group+Sex+race+alq.group,
         nhanes.design, na.rm=TRUE)

##### 2.2 连续性变量 Summary ##### 
# median
svymean(~INDFMPIR+RIDAGEYR+BMXBMI+BMXWAIST+
          total.calories+total.dr.LZ+
          CFDCST1+CFDCST2+CFDCST3+CERAD.total+CFDCSR+CFDAST+CFDDS, 
        nhanes.design, na.rm=TRUE)

svyquantile(~INDFMPIR+RIDAGEYR+BMXBMI+BMXWAIST+
              total.calories+total.dr.LZ+
              CFDCST1+CFDCST2+CFDCST3+CERAD.total+CFDCSR+CFDAST+CFDDS, 
            nhanes.design, na.rm=TRUE,
            quantiles=c(.25,.5,.75),
            interval.type="mean")

#### 3. 课后作业-01 ####
# part1: 衍生 BMI、Smoking status、教育程度的分类变量，命名为 BMI.group, smoke.group, education.attainment，其中 smoking status 可以参考多个定义方式
# part2: 1）膳食补充剂中摄入 L&Z 的平均值的代码，supplement L and Z intake，命名为 total.ds.LZ
#        2）计算总体 L & Z 的摄入量，命名为 total.dr.ds.LZ

